Definition model selection is “estimating the performance of different models in order to choose the best one” (Hastie, Tibshirani, and Friedman 2009)

What it might include

• Choosing structure (GLM, nonlinear model, generalized additive model, discrete/continuous time model, non-likelihood-based [tree, other machine learning methods])
• What variables to include
• How to parametrize (if necessary)

Process in BIOS 650, 651, 653, 675, etc is usually

1. fit model
2. report results

Process in many data analyses is

1. explore data
2. make preliminary decisions
3. fit initial models
4. fit final model
5. report results
6. assess performance

Key idea of this lecture: steps 1–3 impact steps 6 and 7

• \(X \overset{iid}{\sim} MVN(0, \Sigma(\rho))\); \(\epsilon\overset{iid}{\sim} N(0, 1)\)

• \(Y = X\beta + \sigma\epsilon\)

• length of \(\beta\) is \(p\)

• \(n\) observations from model

• \(R^2 = \mathrm{var}(X\beta)/[\mathrm{var}(X\beta) + \sigma^2]\)

• Compare two approaches:

• No variable selection: fit multiple regression with all variables

• Variable selection: exclude all variables with univariable p-value \(> 0.10\), then fit multiple regression with selected variables

library(tidyverse); library(broom);
set.seed(1);
n_sim <- 2e3;
#subjects / simulations
n_subj <- 100;
#compound symmetric correlation
rho <- 0.10;
#true generating betas
p_null <- 19;
#one non-zero beta;
beta <- c(1, numeric(p_null));
#true_rsq implies value of sigma
true_rsq <- 0.2;
#remove all estimates with univariable p-value exceeding 
p.value_threshold <- 0.10;
source("varselect_sim.R");

No variable selection

Variable selection

No variable selection

full_model_coefs %>% 
  filter(term != "x.1") %>%
  group_by(sim_id) %>%
  summarize(num_fd = sum(p.value < 0.05)) %>%
  summarize(mean_num_fd = mean(num_fd), 
            median_num_fd = median(num_fd));

## # A tibble: 1 × 2
##   mean_num_fd median_num_fd
##         <dbl>         <dbl>
## 1       0.904             1

Variable selection

selected_model_coefs %>% 
  filter(term != "x.1") %>%
  group_by(sim_id) %>%
  summarize(num_fd = sum(p.value < 0.05)) %>%
  summarize(mean_num_fd = mean(num_fd), 
            median_num_fd = median(num_fd));

## # A tibble: 1 × 2
##   mean_num_fd median_num_fd
##         <dbl>         <dbl>
## 1       0.518             0

No variable selection

Variable selection

• No variable selection approach unbiasedly estimated regression coefficients but may be overfitting data

• Variable selection approach biasedly estimated regression coefficients but with lower variance and fewer false positives

Some consequences:

• biased coefficient estimates

• non-nominal significance

• false positive findings

• false certainty

• tendency to overstate proportion of variance explained

Variable selection entails a `bias-variance’ trade-off

1. Scientific context should be incorporated wherever possible:

• what is my primary question?

• do prior models exist? how can I incorporate that information?

• are there known constraints? e.g. design features, monotonicity, positivity

1. Data can be used to suggest models

• When used for hypothesis generating, called `exploratory data analysis’ (Tukey 1980)

• SAS / R / STATA have automatic, data-dependent selection routines built in

• When conducted uncritically, ‘data dredging’ or ‘data torture’ are better descriptors

This lecture:

• AIC

• forward / backward selection

Future lecture(s):

• Cross-validation

• LASSO

Variable selection

• Ideal setting: training data to fit model and independent validation data to select model.

• Suppose we have independent observations of outcome given same covariates: \(Y_\mathrm{new} = X\beta + \sigma\epsilon_\mathrm{new}\)

• Akaike (1973) asymptotically linked expected log-likelihood of \(Y_\mathrm{new}\) to observed likelihood of \(Y\)

• Let \(\{\hat\beta_0, \hat\beta,\hat\sigma\}\) be estimated parameters after model building. Then

\[ \begin{aligned} \mathrm{LogLik}_\mathrm{new} &\approx \mathrm{Loglik} - \mathrm{Model~Size}\\ \Rightarrow f(Y_\mathrm{new}|\hat\beta_0, \hat\beta,\hat\sigma) &\approx f(Y|\hat\beta_0, \hat\beta,\hat\sigma) - \#\{\hat\beta\neq0\} - 2 \end{aligned} \]

No variable selection

(small dot indicates mean value)

Variable selection

(small dot indicates mean value)

• AIC defined by \(-2\times\mathrm{Loglik} + 2\times\mathrm{Model~Size}\) (-2 multiplier being for “historical reasons”)

It is…

• metric useful for comparing relative distance of various fitted models to true generating model (smaller is better)

• approximation of expected log-likelihood of your fitted model when sample size is large relative to number of predictors

• more than intuition that subtracting penalty from model fit is a reasonable thing to do

• AIC defined by \(-2\times\mathrm{Loglik} + 2\times\mathrm{Model~Size}\) (-2 multiplier being for “historical reasons”)

It is not…

• absolute measure of model fit

• good approximation of its estimand when sample size is small relative to number of predictors

• useful to compare models fit to non-identical sets of observations

• useful to compare models for different transformations of outcome

• In general, AIC requires likelihoods given ML-estimated parameters, not REML-estimated (see Section 4.7.2, Gałecki and Burzykowski (2013))

• Possible to compare models using ML-estimates, then report REML-estimates from selected model

• See also https://tinyurl.com/yagpx4bv (StackExchange:CrossValidated)

• Hurvich and Tsai (1989) proposed a faster-converging “corrected” AIC with an adjusted estimate of model size:

\[ \begin{aligned} \mathrm{AIC}_C = -2\times\mathrm{Loglik} + 2\times\mathrm{Model~Size}\times\left(\dfrac{n}{n - \mathrm{Model~Size} -1}\right) \end{aligned} \]

• To account for \(p \gg n\), Boonstra, Mukherjee, and Taylor (2015) suggest instead positive part of correction: \(\left(\dfrac{n}{n - \mathrm{Model~Size} -1}\right)_+\)

• Phil’s opinion: always use \(\mathrm{AIC}_C\)

No variable selection

(small dots indicate mean values)

Variable selection

(small dots indicate mean values)

• When considering what variables from a list to include, AIC / AIC\(_C\) are means of comparing them

• But cannot compare all \(2^p\) possible models when \(p\) exceeds 20–30.

This grief has a gravity, it pulls me down But a tiny voice whispers in my mind You are lost, hope is gone But you must go on… Just do the next right thing Take a step, step again It is all that I can to do The next right thing

This grief has a gravity, it pulls me down But a tiny voice whispers in my mind You are lost, hope is gone But you must go on… Just do the next right thing Take a step, step again It is all that I can to do The next right thing

Anna (Frozen 2)

• Forward selection starts with intercept, adds variables sequentially that “improve” model by greatest amount

• Backward selection starts with full model, subtracts variables presence of which “worsen” model by greatest amount

MASS::stepAIC() function does automatic forward/backward selection using AIC as measure of model fit

R functions AIC() and extractAIC() (which stepAIC() uses) employ different constants. From extractAIC help page:

For linear models with unknown scale (i.e., for lm and aov), -2 log L is computed from the deviance and uses a different additive constant to logLik and hence AIC.

For those interested, the difference between AIC() and extractAIC() for linear models is \(n(\log(2\pi) + 1)\)

Digitized image of fine needle aspirate (FNA) of breast mass from 569 patients

Street, Wolberg, and Mangasarian (1993) Mangasarian, Street, and Wolberg (1995)

Outcome Clinical diagnosis (malignant or benign)

Predictors

1. radius (mean of distances from center to points on the perimeter)
2. texture (standard deviation of gray-scale values)
3. perimeter
4. area
5. smoothness (local variation in radius lengths)
6. compactness (perimeter^2 / area - 1.0)
7. concavity (severity of concave portions of the contour)
8. concave points (number of concave portions of the contour)
9. symmetry
10. fractal dimension (“coastline approximation” - 1)

Data include mean, standard error, and worst measurements.

Available at https://archive.ics.uci.edu/ml/datasets/breast+cancer+wisconsin+(diagnostic)

breast_dx <-
  read_csv("bdiag.csv") %>%
  # Translate M/D into 1/0
  mutate(malignant = 1 * (diagnosis == "M")) %>% 
  # Drop errant space in 'concave points_mean' variable name 
  rename_with(~str_replace(string = ., pattern = " ", replacement = "")) %>%
  # Focus only on worst measurements
  select(malignant, 
         #contains("_mean"), 
         #contains("_se"), 
         contains("_worst")) 

Plan to fit a logistic regression model to predict probability of being malignant/benign using set of worst measurements

ggcorrplot(cor(select(breast_dx, -malignant)), hc.order = T)

1. Sample equal sized dataset with replacement from the original dataset
2. Run forward selection process bootstrap dataset, separately for each criterion
3. Record selected model
4. Repeat steps 1–3 \(B=2000\) times

Total of 142 distinct variable combinations selected across the B = 2000 bootstrap datasets

AIC, AIC\(_C\) require two ingredients

\(\mathrm{AIC}\), \(\mathrm{AIC}_C\) useful model selection technique depending upon being able to do two things:

• calculate likelihood

• count parameters (not always trivial)

When one or both of these is not calculable, must use other options, e.g. cross-validation / split-sampling

There are better options than stepwise selection

• Least absolute shrinkage and selection operator (LASSO)

• Elastic net

Be familiar with your data / interrogate your model

Model selection is not model assessment

Focus of talk has been on model selection (“estimating the performance of different models in order to choose the best one”), distinct from model assessment (“having chosen a final model, estimating its prediction error…on new data”, Hastie, Tibshirani, and Friedman (2009))

Standard errors of selected model do not account for selection process. Ideal assessment requires additional (third) data set: testing data (future lecture (?))